Nonequilibrium molecular dynamics of complex fluids near the gel point 
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We have carried out nonequilibrium molecular dynamics simulations of a system of crosslinked 
particles under shear flow conditions. As the fraction of crosslinks p is increased the system ap- 
proaches a gel point at which the shear viscosity n and the first and second normal stress coefficients 
and ^2 diverge. All three quantities seem to diverge with a power law form: n ~ e~ s , 9i t 3 ~ e~ x 
where e = p c — p and s w 0.7 and A tn 3.15 in three dimensions and s w 2.0 and A m 6.0 in two 
dimensions. 

PACS numbers: 61.43.Hv,66.20.+d,83.10.Rs,83.60.Bc 
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I. INTRODUCTION 

There have been many studies, both experimental and 
theoretical, of the rheological properties of complex fluids 
in the vicinity of the gel point Q] . There is general agree- 
ment in the literature that the transition from the sol to 
the gel phase, at least in the case of chemical gels, is con- 
tinuous and accompanied by the divergence of the shear 
viscosity 77 ~ (p c — p)~ s where p characterizes the de- 
gree of crosslinking or condensation and p c is the critical 
point. However, the experimental and theoretical values 
of the exponent s are rather widely distributed and there 
is continuing debate concerning the existence of a sin- 
gle universality class for these transitions. We mention, 
in passing, that we are excluding the case of vulcaniza- 
tion from this discussion. There is substantial evidence, 
in the form of a Ginzburg criterion 0, that the vulcan- 
ization transition which involves the crosslinking of very 
long chains is for practical purposes mean-field-like. 

There has been much less theoretical and experimen- 
tal work on the normal stress differences. In a non- 
Newtonian fluid under shear flow, e.g., in a Couette 
geometry with the flow in the ^-direction and velocity 
gradient in the ^-direction, the first and second normal 
stress differences iVi = a xx — a zz and N2 = o zz — a yy 
are both not zero and, in the low shear-rate limit, pro- 
portional to the square of the shear rate 7. It is con- 
ventional to define the normal stress coefficients ^2 
through Ni = 1 Fi7 2 . One of the purposes of this work is 
to investigate the critical behavior of the normal stress 
coefficients. These are most easily calculated directly, 
i.e., by shearing the computational cell and measuring 
the normal stresses via a virial formula. While there ex- 
ists a Green-Kubo formula Q that yields the zero shear- 
rate limit of the integral involved presents significant 
computational problems close to the gel point: the decay 
of the relevant correlation function becomes extremely 
slow and truncation errors become unmanageable. For 
this reason, we have decided to utilize nonequilibrium 
molecular dynamics (NEMD) to calculate both the shear 



viscosity and the normal stress differences. We use the 
same model for which we have previously calculated the 
shear viscosity Q using the Green-Kubo formalism. Our 
NEMD results for the critical exponent s are consistent 
with these earlier results. The NEMD simulations indi- 
cate that the first normal stress coefficient diverges at 
the critical point ^1(7 = 0) ~ (p c — p)~ x with A ~ 3.15, 
i.e., much more strongly than does the shear viscosity. 
The errors in ^2(7 = 0) are much larger than those in 
'Fi and an independent determination of an exponent 
for its divergence is not feasible. However, the data are 
consistent with the conjecture that both normal stress 
coefficients diverge in the same way. We note that the 
conclusion that 'Fi diverges more strongly than 77 has also 
been arrived at by Broderix et aL[3j] in the context of a 
Rouse-type model. However, their exponent A ~ 4.9 is 
significantly larger than ours. 

The structure of this article is as follows. In section 
[H]we briefly describe the model that we have used and 
the computational techniques. Section IlIII contains our 
results for the two-dimensional case and the results in 
three dimensions are presented in Section IIVI We con- 
clude in Section W\ with a brief discussion. 



II. MODEL AND COMPUTATIONAL 
METHODS 

Our model of the sol phase is identical to that in Q| . 
All particles interact through the soft sphere potential 

V{ r ij) = e ( f7 / r y) 36 f° r r ij — l-^c an d, in the three- 
dimensional calculations, we have used a single volume 
fraction $ = irNa 3 /6V = 0.4 which is well below the 
liquid-solid coexistence density. All calculations were 
carried out at a temperature ksT/e — 1. In the ab- 
sence of crosslinks, this system is a simple liquid that 
has been well characterized |g. We initially placed the 
particles on the vertices of a simple cubic lattice and in- 
stantaneously and randomly introduced a fraction p of 
nearest-neighbor bonds. We used the bonding poten- 
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tial Vb(rij) — k(rij — ro) 2 where k = 5e/<7 2 and where 
ro = (7r$) 1 / 3 /cr so that there was no internal mechani- 
cal strain. This method of crosslinking ensures that the 
cluster size distribution is that of percolation in three di- 
mensions and that a gel forms (in the thermodynamic 
limit N — > oo) at p c w 0.2488. Once the particles 
had been crosslinked they were free to move through- 
out the three-dimensional computational box. They were 
initially thermalized with periodic boundary conditions. 
Once equilibrium had been attained, the computational 
box was sheared at a rate 7 = dv x (z) jdz and the bound- 
ary conditions were changed to the Lees-Edwards bound- 
ary conditions 0. The system was then reequilibrated 
using the so-called SLLOD algorithm Q subject to the 
constraint that the kinetic energy in the frame following 
the overall flow of the particles remain constant. This 
kinetic energy is proportional to the square of the "pe- 
culiar velocity" (v x — A fz,v y ^v z ). Once a steady drift had 
been established, the diagonal elements of the stress ten- 
sor as well as a xz were calculated from the appropriate 
virial formula. Calculations were performed for systems 
of N = L 3 particles with L = 10, 15 and 20 over the 
entire range < p < p c and for shear r ates ge nerally in 
the range .005 < 7 < 0.1 in units of ^JTj ma 2 although 
for a few cases larger shear rates were also imposed. It 
should be noted that if atomistic values of e, m and a 
are used, this range of shear rates corresponds to values 
of order 10 12 s~ 1 , i.e., enormously large rates compared 
to experimental values. Even if one takes the point of 
view that the particles represent a colloidal suspension, 
the minimum shear rate is still of the order of 10 3 s _1 . It 
was necessary to use such large shear rates in order to ob- 
tain reasonably well converged estimates for the normal 
stress coefficients, especially close to the gel point. For 
this range of shear rates, the equations of motion could 
be stably integrated with a time step St = 0.005 y / mcr 2 /e 
but for larger values of 7 the time step had to be de- 
creased. 



The values of the shear viscosity, and to an even larger 
extent the normal stress coefficients, varied considerably 
for different realizations of the crosslinks. Therefore, we 
averaged the results over several thousand realizations 
even for values of p as small as 0.1. 



The same repulsive pair potential was used in the two- 
dimensional case. Here the particles were initially placed 
on a triangular lattice and instantaneously crosslinked 
as in three dimensions. This sets the gel point at p c — 
2 sin(vr/18) « 0.347296. The lattice constant of the initial 
configuration was 1.2, and so these simulations were done 
at a number density n « 0.8. The spring constant used 
was k = 40e/cr 2 and the value of ro = 1.2 was chosen to 
eliminate mechanical strain due to the crosslinks at zero 
temperature, as was done in three dimensions. 
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FIG. 1: The viscosity as a function of p for a crosslinked fluid 
in two dimensions. The straight line shows a power law, with 
exponent s — 2. The inset shows the shear rate dependent 
viscosity (open symbols) as compared to the zero shear rate 
values (filled symbols). The curves show fits of the finite shear 
rate data to the Lorentzian form discussed in the text. 



III. RESULTS IN TWO DIMENSIONS 

The final results for the zero shear rate viscosity ex- 
trapolated from nonequilibrium molecular dynamics sim- 
ulations are shown in Fig.^ These data are from simula- 
tions of systems of size 32 x 32 = 1024 particles. No sig- 
nificant deviation from a scaling form is visible for p > 
at this system size for the range of p studied here. The 
power law divergence of the viscosity thus directly gives 
an estimate for the exponent s = 2, without a finite size 
scaling analysis. 

The inset to Fig. ^ shows the shear rate dependence 
of the viscosity as well as the zero shear rate values from 
the Green-Kubo formula, for an uncrosslinkcd and for a 
lightly crosslinked sample. The extrapolation of the fi- 
nite shear rate values to zero shear-rates seems to be con- 
sistent with the Green-Kubo values. This fluid exhibits 
shear thinning, as is seen in some experiments on complex 
fluids. The determination of a zero-shear-rate viscosity 
requires fitting to some functional form for the shear- 
rate dependence of rj. However, this value is not sensi- 
tive to the form chosen. Many different functional forms 
have been suggested for the shear-rate dependence of the 
viscosity, mostly suggested as phenomenological fitting 
functions to experimental data p|. The value for 77 in 
the main figure was estimated from a fit to a Lorentzian 
plus a constant term, as suggested in 0, while the dif- 
ference between different fits was used as an estimate of 
the error. The Lorentzian form has the advantage that 
it is automatically symmetric in 7 and is analytic near 
7 = 0. 

There is only one normal stress difference, N\ = a xx — 
(Tyy, in a two dimensional system. We have measured 
the associated coefficient ^1 close to p c . This quantity is 
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FIG. 2: The normal stress coefficient $ in two dimensions 
close to p c . The straight line is a power law with exponent A = 
6. The inset shows the linear extrapolation used to determine 
the 7 = value of 



expected to diverge as a power law 'J'i ~ {p c —p)~ x as the 
gel transition is approached; the results of our simulation 
are shown in Fig. [21 We estimate an exponent A « 6 from 
our data. 

The estimation of the normal stress coefficient was 
more difficult than for the viscosity. The division by -y 2 
produces large statistical errors for small 7, thus making 
it difficult to determine the functional form of the shear- 
rate dependence. An example is shown in the inset to 
Fig. |21 Only a simple linear fit was used to estimate the 
zero-shear-rate normal stress coefficient. 
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FIG. 3: Dependence of the shear viscosity r/(p, 7) on the shear 
rate 7 for L = 10 and a selection of crosslink densities. 
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IV. RESULTS IN THREE DIMENSIONS 

In the three-dimensional case, we have results for sys- 
tems of N = L 3 particles with L = 10, 15, and 20. 
We first display, in Fig. |3 the shear-rate dependence 
of the viscosity 77(7) for L = 10 and five different val- 
ues of the degree of crosslinking ranging from the sim- 
ple fluid case, p = 0, to a system close to the gel 
point (p = 0.22). All systems show evidence of shear 
thinning, with this feature becoming much more promi- 
nent and setting in at lower shear rates as the gel point 
is approached. If we rescale the shear viscosity using 
the form 77(^,7) = a(p c — p) s f){p, 7) with ,s = 0.7 and 
7 = b(p c — p)~ ZA { where a and b are constants, we can 
achieve a respectable collapse of the data, as seen in 
Fig. 0] Since the data are noisy, we have not made a 
serious effort to optimize this collapse. Nevertheless, the 
"dynamical" exponent z w 2.35. The exponent s = 0.7 
used to rescale the viscosity is our best estimate of the ex- 
ponent that governs the divergence of the zero shear-rate 
viscosity 77(^,7 = 0) ~ (p c -p)~ s - 

We note that dynamical scaling yields a connection be- 
tween the exponent z that controls the divergence of the 



FIG. 4: The data of figure 1 plotted as function of the scaled 
shear shear rate 7 = (p c —p)~ ZA f with 2 = 2.35. The viscosity 
has been multiplied by (p c — p) s with s — 0.7 to remove the 
divergence at p c . 

longest relaxation time at the gel point and the expo- 
nents s and t, i.e., z = s + 1 1]. Here t is the exponent 
that describes how the shear modulus vanishes as the gel 
point is approached from the solid side: /1 ~ {p — Pc) 1 - 
For this model, we have determined ^3] that t w 2.0. 
This yields the prediction z w 2.7, a value not too far 
from the value used to rescale the shear rate. 

The zero-shear-rate viscosity is shown in Fig. in 
finite-size scaled form, i.e., we plot L~ s l v r\(j>, L) as func- 
tion of (p c — p) p L where v is the correlation length 
exponent of the three dimensional percolation problem 
v k 0.88. As mentioned above, the value of the exponent 
s = 0.7 is consistent with our previous result obtained 
from a Green-Kubo calculation 

We turn now to the normal stress coefficients. In 
Fig.Elwe plot the first normal stress coefficient ^1 (p, 7) = 
(&XX — o'zz)/'! 2 as function of the shear rate for L = 10 
and three values of p. For larger values of p, vFi increases 
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FIG. 5: Finite size scaling plot of the zero-shear-rate viscosity. 



rapidly as the shear rate is decreased and an estimate of 
the zero shear-rate value 7) is problematical. We 

have fit the data points to a second order polynomial 
^i(p>7) — ^i(p> 0) + (17 + fr-y 2 and obtained our estimate 
of the zero shear-rate value in this way. A fit to an expo- 
nential decay works equally well and produces estimates 
of \&i(p, 0) that differ by no more than 3% from those 
shown here. Similarly, the Lorentzian-plus-constant fit 
that was used to fit the viscosity in two dimensions also 
provides a reasonable fit to the data as long as there are 
enough values of the shear rate (more than five). The 
conclusions presented below are insensitive to the method 
of extrapolation. 

The unsealed data for *f?i(p, 0) are plotted in Fig. for 
L = 10, 15 and 20 as function of p c — p along with a line 
representing the function a(p c — p)~ x with A = 3.15 that 
captures the form of the data outside the critical region 
quite well. Closer to the critical point p c « 0.2488, the 
usual finite-size effects that appear when the geometric 
correlation length approaches the system size, L, are 
evident. These finite-size effects can be hidden by plot- 
ting L~ x '"tyi(p, L) as a function of the scaled variable 
e = L/t; oc L(p c — p)" where v « 0.88 is the correla- 
tion length exponent. This is done in Fig. [8] and the 
data do collapse reasonably well to a universal curve. A 
useful consistency check on this procedure is available 
far from the critical point: if the exponents A and v are 
correctly determined then the scaled normal stress coeffi- 
cient should approach the power-law form £~ x f u at large 
e. This line is also displayed in Fig. [S] and the data are 
consistent with the expected behavior. 

We have also calculated the second normal stress coef- 
ficient ^2(p, L). This coefficient is much smaller in mag- 
nitude than \T/i and negative, at least for our range of 
shear rates. The sample-to-sample fluctuations in N2 are 
relatively much larger than those in N\ and the data for 
SS?2 are therefore much more noisy. Indeed, they are so 
poorly converged that a determination of a critical expo- 
nent from that dataset is not supportable. We have there- 
fore carried out a rescaling of \&2 under the assumption 
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FIG. 6: First normal stress coefficient for L = 10 and p = .15, 
.2 and .21. The curves are fits of the data to a second order 
polynomial. Note the rapid increase of $1 (p, 7 = 0) as p — > p c . 




FIG. 7: The first normal stress coefficient ^i(p, L) plotted as 
a function of p c — p. The solid line corresponds to a power 
law divergence of the form ^i(p, 00) ~ (p c — p) -3 ' 15 . 
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FIG. 8: The first normal stress coefficient multiplied by L~ x ^ 
with A = 3.15 plotted as a function of e = L(p c — p)" ■ The 
straight line correspond to the expected large-e form e -A//l 7 
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FIG. 9: Finite size scaling plot of the second normal stress 
coefficient normal stress coefficient using the same value of A 
as in Fig. |H] 

that it diverges in the same way as $i, i.e., controlled 
by the same critical exponent A w 3.15. The results of 
this rescaling are shown in Fig. and we see that the 
data for L = 15 and L = 20 that are somewhat better 
converged than the data for L = 10 seem to support such 
an assumption. 

The authors of reference || have proposed a scaling 
relation for the exponent A: A = z + s that seems to be 
rigorous for the Rouse model that they have used. In 
our case, using the two estimates of z referred to above, 
namely z = 2.35 and z = 2.7, we obtain A = 3.05 and 
A = 3.4 which bracket the measured value. It should 
however be noted that in the Rouse model of j3j, the 
second normal stress coefficient ^2 = for all values of 
P- 



ity and the normal stress coefficient in a model gel as 
the gel transition is approached. Our results for the di- 
vergence of the zero shear-rate viscosity are consistent 
with our previous calculation using a Green-Kubo for- 
mula to extract the viscosity from an equilibrium simu- 
lation. In addition, this model exhibits shear thinning as 
the shear rate is increased, as is observed in experiments 
on gelling systems. We found that the exponent govern- 
ing the divergence of the shear viscosity to be s = 0.7 
in three dimensions. This value is consistent with some 
experiments 0, 0] , as well as several analytical calcu- 
lations We also find s — 2 in two dimensions. 

We have also presented evidence that, in this model, 
the shear-rate-dependent viscosity can be rescaled onto a 
single universal curve. This indicates that the physics of 
the shear-thinning that we observe for all p is the same 
close to the gel point as it is in the simple liquid (p = 0). 

The divergence of the normal stress close to the gel 
transition has not been measured in a experiment. As 
suggested in , it should be possible to observe the very 
strong divergence of this quantity experimentally. Mea- 
suring both the viscosity and the normal stress close to 
the gel point would give the ratio of two dynamical ex- 
ponents without determining the critical point, which is 
often difficult to determine accurately in an experiment. 
The ratio of s to A would then provide a dynamical expo- 
nent which would characterize the universality class of a 
given material. Comparing this value to the values pre- 
dicted by different models could then give some insight 
into which features of a microscopic model are important 
to the dynamical properties of an incipient gel. 
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